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INTRODUCTION 


This  note  is  concerned  with  a problem  of  numerical  evaluation 
that  arises  in  the  analysis  of  oscillating  lifting  surfaces;  its  purpose  is 
to  indicate  a method  of  obtaining  accurate  solutions  where  no  such  method 
seemed  to  exist,  and  to  illustrate  the  magnitude  and  the  nature  of  the 
errors  that  arise  from  current  approximate  methods.  Completeness  is 
not  intended;  that  is,  while  given  problems  may  involve  various 
complications,  depending  on  the  model  at  hand,  those  complications  that 
are  not  relevant  for  the  intended  demonstration  of  a basically  different 
approach  will  not  here  be  considered.  Rather,  we  assume  the  simplest 
model,  the  planar  wing  in  incompressible  flow.  (The  planform  of  the 
wing  is  not  of  concern.  ) 

Also,  while  the  new  approach  may  be  generalized,  we  confine  our 
discussion  to  a specific  finite  element  method  of  solution,  namely,  to  the 
Falkner  method.  In  finite  element  methods,  generally,  the  surface  is 
treated  as  a composite  of  finite  elements  Er;  one  requires  the  aerodyna- 
mic  influence  coefficient,  which  measure  the  downwash  due  to  unit 

(oscillating)  lift  on  any  "sending"  element  Eg,  on  pre-determined  colloca- 
tion points  (xn,  yn)  on  the  "receiving"  elements  En  (which  include  Eg). 

In  the  Falkner  "horseshoe  vortex"  method  (also  called  "doublet- 
lattice"  method,  see,  e.g.,  Ref.  1)  the  sending  element  Eg  is  modeled 
by  a horseshoe  vortex  the  bound  part  of  which  (i.  e.  , the  part  that  pro- 
duces the  lift)  spans  the  l/4-chord  line  of  Eg.  For  this  model,  we  will 
first  give  an  indirect  (and  tentative)  indication  of  the  likely  distribution  of 


2 


* 

W \ . 


the  errors  which  arise  from  a specific  approximative  method,  here  to 


1 


be  called  the  "parabolic  approximation"  . We  will  then  show  how  the 


wg  n can  be  calculated  to  arbitrary  accuracy,  and  will  demonstrate 


numerically,  for  the  two  most  critical  w , the  errors  that  exist  at  two 

s i n 


different  levels  of  approximation.  Finally,  we  will  briefly  discuss  the 
status  of  the  overall  problem  and  its  remaining  difficulties. 


FORMULATION 

We  will  use  exclusively  the  configuration  of  elements  which  is 


shown  in  Fig.  1.  The  are  equal  squares.  The  x-axis  is  parallel  to 


the  undisturbed  flow  speed,  V,  and  is  centered  on  the  sending  element  Ec 


The  bound  vortex  lies  on  the  y-axis.  Both  coordinates,  x and  y,  are  re- 


ferred to  the  length  i ; 2i  is  the  side  length  of  each  En-  The  same  length 


f we  use  to  define  the  reduced  frequency 


k = Wf/V  (1) 

By  unit  lift  (above)  we  imply  that  the  lift  coefficient  of  Eg  is 

C = 1 • exp  (i«t)  (2) 

J-j , s 


The  influence  coefficient  w measures  the  dimensional  downwash  W 

s,  n s 


at  the  collocation  point  (xn>  yn)  as  follows: 


it  W (x  , y )=w  V-  exp  (Lot) 
s n ■'n  s,  n c 


(3) 


Though  we  will  require  only  individual  influence  coefficients 

w , these  are  really  individual  values  of  a continuous  (in 
s,  n 

general)  function  w (x,y;k),  and  in  our  analysis  it  is  more  convenient  to 


: 


■ f 


. j 

i 


i 


i 


j 


use  this  functional  form.  With  the  definitions  introduced  above,  this 

influence  function  is  given  by  an  integral  of  the  form 

+1 

ws(x,y;k)  = J K(x,|y-n  |,k)  (4) 

All  the  terms  in  Eq.  (4)  are  non-dimensional.  The  detailed  form  of  the 

kernel  K we  will  discuss  later. 

Whenever  y=0  (i.  e.  , when  En  is  centered  on  the  x-axis)  the 

integral  in  Eq.  (4)  is  a Hadamard  integral,  singular  to  second  order.  On 

the  other  hand,  the  kernel  K is  given  numerically  and  to  limited  accuracy. 

For  this  reason,  presumably,  the  scatter  between  results  obtained  by 

2 

various  current  methods  is  around  10%.  This  is  not  satisfactory;  our 
goal  is  to  obtain  accurate,  reliable  results. 

THE  PARABOLIC  APPROXIMATION 

In  this  section  we  want  to  illustrate,  though  in  an  indirect  manner 
(by  means  of  readily  available  exact  results)  the  working  of  a current 
approximate  method.  Our  purpose  is  to  obtain  a tentative  indication  as  to 
which  parameter  range  should  be  of  foremost  concern. 

In  the  steady  limit  (k=0);Eq.  (4)  can  readily  be  resolved  analyti- 
cally, since  here  the  kernel  takes  the  simple  form 

= 1 + x/ (x^+  (y-r|)^  )2  (5) 


The  suffix  (s)  in  K ^ denotes  "steady". 


r 

In  the  method  proposed  in  Ref.  1,  the  unsteady  part  K'u'  is 
defined  by 

K = K(s)+  K (u)  (6) 

This  rather  complicated  part  is  approximated  by  the  parabola  through 

its  values  for  r|  = -1,  *1=0  andq  = l.  The  integration  over  this  parabola  is 

performed  analytically,  and  the  "parabolic  approximation"  to  the  unsteady 

part  w ^ (x,  y)  of  the  influence  function  is  thus  obtained.  By  adding  the 

known  exact  contribution  of  K'  ' , one  obtains  an  approximate  value  for 

the  complete  influence  function  w (x,  y;k). 

s 

( g ) 

We  apply  here  the  parabolic  approximation  to  K (rather  than  to 
K^)  and  can  then  compare  this  approximation  with  the  known  exact  re- 
sult. For  example,  for  y=0  the  approximation  is 

2w jS^(x,  0-,0)  = | - x/(HxV  (x<°)  (7) 

and  the  exact  result  is 


Numerical  values  are  listed  in  Fig.  1 (where,  as  is  characteristic  for 
the  Falkner  method,  the  central  3/4 -chord  point  of  each  En  is  used  as 
its  collocation  point).  The  correct  result,  Eq.  (7a),  is  written  above 
each  collocation  point,  the  error  Awg  approximation  Eq.  7 is 

written  underneath.  A bar  (-)  underneath  indicates  that  the  (absolute) 


error  is  <0. 0001 . 
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The  errors  listed  in  Fig.  1 (antisymmetric  with  respect  to 
the  y-axis)  converge  to  zero  rapidly  away  from  Eg  --  much  more 
rapidly  than  the  wg  themselves  converge  toward  their  respective  limits. 
Sizable  errors  arise  on  and  near  Eg;  the  largest  errors  arise  along  the 
x-axis  (where  the  integral  in  Eq.  (4)  is  singular).  Along  this  axis, 
exact  influence  function  and  approximation  are  shown  in  Fig.  2,  where  it 
is  seen  that  the  error  becomes  infinite  as  the  bound  vortex  is  approached. 
In  Fig.  1,  the  largest  percentage  error  is  29%  (at  the  collocation  point 

(-1,0))  and  is  thus  rather  large. 

Away  from  the  x-axis,  the  errors  listed  in  Fig.  1 are  much 
smaller,  an  important  reason  being  that  outside  the  bound  vortex  the 
error  at  x=0  is  zero  rather  than  infinite.  The  largest  errors  along  the 
verticals  y = i 2 are  A w ^ = 'tO.  0009  at  x - _0.  5. 

AN  ACCURATE  METHOD 


a.  General 

In  the  case  of  unsteady  flow,  the  kernel  K is  given  by  the  integral 
K = J exp(-ik(x+|y-p|  u))  TITT  ^ 

-yiUi  d*'*2)  ' 

To  evaluate  Eq.  (8)  directly  means  that  one  has  to  deal  with  incomplete 
cylinder  functions.  The  evaluation  method  proposed  in  Appendix  A 
of  Ref.  1 is  not  particularly  accurate.  There  are  thus  two  distinct 
sources  of  possible  error:  the  evaluation  of  K,  and  the  parabolic  approxi- 

*A  new,  slightly  simpler  but  much  more  accurate  method  is  pre- 
sented in  Ref.  3. 


6 


mation.  One  would  tend  to  expect  that  the  second  source  would  produce 

the  larger  errors  and,  therefore,  would  tend  to  expect  that  the  largest 

errors  will  again  arise  at  the  collocation  points  ( - 1,  0)  as  in  Fig.  1.  On 

this  basis  and  in  order  to  simplify  the  presentation  of  the  essence  of  our 

method,  we  assume  in  the  present  section  that  y = 0.  For  brevity,  then, 

we  write  w (x;k)  for  w (x,  0;k).  (In  fact,  as  we  will  see  later  (Eq.  18), 
s s 

wg(x,  y;k)  can  always  be  expressed  by  wg(x;k). ) 


Equation  (4)  with  Eq.  (8)  inserted  presents  a double  integral. 

The  essence  of  our  approach  is  that,  by  reversing  the  sequence  of  inte- 
gration, we  avoid  the  incomplete  cylinder  functions.  (That  this  reversal 
is  permitted  one  checks  most  readily  in  the  steady  case.  Here  all  inte- 
grations can  be  performed  analytically  in  either  sequence,  and  the  same 
result  is  obtained  either  way.  ) 

Setting  y=0,  ur|  = £,,  and  reversing,  one  obtains 

1 


-2w  (x 
s 


co  1 00 

’k)eikx  = f e-*t(  / -2d\  3/2  W = / f— " J< 
J (/  (on /z  -x  £2(i+n2 


■X)  (9) 


The  problem  of  determining  wg  has  been  reduced  to  that  of  evaluating  the 
single  (complex)  integral  J(-x).  The  real  and  imaginary  parts  of  this 
integral  we  will  denote  as  follows 


J(x)  = JC(x)  - iJS(x) 


(9a) 


the  suffix  c indicating  cos,  the  suffix  s indicating  sin  in  the  integrand. 

The  integral  J is  not  difficult  to  evaluate,  tut  different  methods  are 
required  in  different  ranges  of  the  variable  t,.  To  indicate  an  integral 
that  extends  from  £,  = a to  £,  = b,  we  write  J(a,b),  but  we  usually  omit  the 
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upper  limit  if  this  limit  is  °o.  Thus 


J(a)  = J (a,  oo)  = J(a,  b)  + J(b) 


The  singularity  at  f,  = 0 enters  if  x is  positive  but  not  if  x is 
negative,  and  this  turns  out  to  be  a major  distinction.  In  order  to  avoid 
confusion  of  signs  (and  a clumsy  notation),  we  introduce  p to  mean  (xl  . 
We  deal  first  with  the  case  x ^ - 1 where  there  is  no  singularity. 


b.  The  Case  x<  -1 


Writing  v for  k£,,  we  have 


7 r 00 

(p)=k1 2  / COSV_^  dv^y  jC( 

4 (1  +(k/v)2)  2 v3 *  v 

oo 

-C,  . , ,v-l . 2v  r cosvdv 

with  J (p)  = (-  ) v-lk  / 2v+l 

V 7 v 

kp 


(10a) 


1 • 3 • 5 ■ • (2v  - 1 ) - 

where  v - Z A b - (2v ) : 


(10b) 


The  same  formulas  apply  to  J(p)  if  cos  is  replaced  by  sin.  The  two 

types  of  elementary  integrals  are  given  in  Ref.  4,  formulas  333:7b  and 

6b.  After  rearrangement,  they  take  the  forms 


£Jc(p) 


= Acos(kp)  - kBsin(kp)  +k  CCi(kp) 


(Ha) 


|jS(p)  = Asin(kp)  + kBcos(kp)  + k Csi(kp) 


(lib) 


The  functions  Ci  and  si  = Si  -ir/2  are  the  cosine  and  sine  integrals  (e.g., 
Ref.  5,  formulas  5.  2.  1, 2&5).  These  may  be  calculated  to  very  high 
accuracy  by  means  of  IBM  Subroutine  SICI,  Ref.  6,  p.  370. 
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Insertion  of  Eqs.  (11a,  b)  into  Eq.  (9)  yields 


w (x;k)  = -A+ikB-k  [ Ci(kp)  - is  i(kp)  ] Ce  (x^~l) 


(12) 


If  the  argument  kp  is  large,  Ci  and  si  are  best  expressed  by  the  functions 
f(z)~l/z,  g(z)~l/z^ , see  formulas  5.  2.  8,  9,  34  & 35  of  Ref.  5.  Then 


2 

Wg  (x;k)  = -A  + ikB  -k  [ if(kp)  - g(kp)]C 

No  oscillatory  term  is  left  in  this  asymptotic  form. 

The  functions  A,  B and  C of  Eq.  (11)  are  given  by  series: 


(12a) 


A = £ k2nA 

o n 


2 n. 


1 00  2n  00  A 

Z k Z (-) 

2p2  o V=D  v 


2n 


vB 


B = 2 k B = 2 2 (-) 

o n 2P  o v = o 2 


C = 2 k2nC  = i ^k2n  n 


n,  v 
>v 


(2n  + 2)! 


a _ n+v  (2v) ! 

A = (2v+l)B  ; B = pz — 7—5 — 7—5-.  1 

n,  v n,  v n,  v (2n  + 2v+2)! 


(13a) 


(13b) 


(13c) 


(13d,  e) 


Note  that 


2 A = -1+  (l+p2)Vp 


(13f) 


If  one  sets  k = 0 in  Eq.  (12),  all  terms  on  the  right  disappear  except  -Aq 
and  the  steady  flow  result  Eq.  (7a)  is  thus  recovered. 

Noting  further  that 

A n = B n = C 
n,  0 n,  0 n 


( 1 3g) 


one  sees  that  one  can  simplify  by  splitting  off  the  terms  v = 0, 
setting 
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l 


A = C/p2-A;  B = C/  p - B (13h,  k) 

so  that  Eq,  (12)  becomes,  again  for  x 5 - 1, 

wg(x;k)  = - { l/p2  - ik/p  + k [ Ci(kp)  - isi(kp)]  e^^}  C+A-ikB  (14) 
The  factor  C is  quickly  calculated: 

1 k^  k^  k^ 

C=  4+9^+7840  + 255048  + • • • (14a) 

The  remainder  terms  we  write  as 

00  ® 'y 

A = £ k2n  A (P)  ; B=  D k nBn(p)  (14b,  C) 

Sample  numerical  values  are  given  in  Table  1;  using  these,  one  can 
quickly  find  A and  B if  p = 1 , 2 or  >3.  For  p -values  in  between'",  only 
Aq  and  Bq  may  be  required  if  k is  small  enough  and/  or  if  p is  large 
enough  . Of  these,  Aq  can  be  calculated  from  its  analytical  form,  com- 
pare Eq.  ( 1 3 f ) 


Aq  = [1  + 1/2  p2  - (1  + 1/p2) M/2  (14d) 

On  the  other  hand,  the  series  Bq  converges  rather  slowly  if  p is  close  to  1. 
Fortunately,  one  has  an  easy  and  quick  alternative  in  all  the  cases  when 
evaluation  of  the  double  series  becomes  inconvenient;  returning  to  Eq.  (9), 
one  calculates  the  part-integral  J(p,  b)  (with  either  b=2  or  b = 3)  by  numerical 
integration,  and  determines  the  remaining  integral  J(b)  using  Table  1. 


5? 

In  Fig.  1,  x is  always  an  odd  interger.  Other  values  of  x do  arise, 
however,  e.  g.  in  compressible  flow,  orwheny  / 0,  see  Eq.  (18)  below. 
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c.  The  Case  x > + 1 

If  x ?+  1,  then  x = p,  and  from  the  symmetry  properties  of  the 
circular  functions  follows  that 


Jc(-p)  = 2JC(0)  - JC(p)  ; JS(-p)  = JS(p)  (15a.  b) 

Therefore,  Eq.  (9)  takes  the  form 

w (x;k)eikx  = |jC(x)  - JC(0)  JS(x)  (x>+l)  (16) 

s ^ 


An  alternate  form  of  Eq,  (16)  is  sometimes  useful: 


w (x;k)  = -Rew  (-x;k)  + ilmw  (-x;k)  - JC(0)e 
s s s 


(16a) 


Of  the  terms  on  the  right  of  Eq.  (16),  only  the  part-integral  J (0,  1)  has 

not  already  been  disucssed  in  the  preceding.  For  this  integral  we  have 

1 ^ 

coskg,dt,  y->v2vT 

(17) 


-JC(0,  1)  * - f -- 


The  elementary  integrals 
1 


I = 
v 


(-)v  + 1 f t2v-2dt 


(2v)! 


,/L 

o <* 


d+r) 


2,  i 


[ 17a) 


are  standard  integrals  (e.g.,  Ref.  4,  section  234).  Their  numerical 

values  converge  to  zero  rapidly  as  v increases;  see  Table  2. 

c 

When  x is  large  (and  positive),  the  integral  J (0)  provides  by  far 
the  largest  contribution  to  wg(x;k),  and  this  contribution  is  oscillatory. 
Herein  lies  a pronounced  distinction  between  the  ranges  of  large  positive  x 
and  of  large  negative  x,  Eq.  (12a). 
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(d)  The  Case  x < 1 

c s 

In  this  case  one  requires  the  part-integrals  J (p,  1)  and  J (p,  1). 
Series  presentations  corresponding  to  Eqs.  (17),  ( 17a)  can  be  used 
but  will  be  less  convenient  then  direct  numerical  integration  (except  if 
k is  quite  small).  Equations  (16)  and  (16a)  remain  applicable. 

NUMERICAL  EXAMPLES 

The  largest  errors  from  applying  the  parabolic  approximation  to 
the  steady  flow  problem  occurred  (see  Fig.  1)  at  the  collocation  points 
closest  to  the  bound  vortex,  that  is,  at  the  points  (x,  y)  = (t  1,  0).  For 
these  points  we  have  calculated  (for  the  range  0 ^ k ^1)  the  following  sets 
of  results: 

(a)  the  accurate  influence  coefficients  w ; 

s 

(b)  the  wg  calculated  by  applying  the  parabolic  approximation 

to  only  the  unsteady  part  of  the  kernel,  using  accurate 

values  for  K^; 

(c)  the  corresponding  parabolic  approximations  calculated  using 

the  less  accurate  values  for  that  one  obtains  using  the 

method  proposed  in  Appendix  A of  Ref.  1. 

In  Figs.  3 to  5 some  of  these  results  are  presented  and  are  denoted  by 
"correct",  "parabolic  approximation",  and  "Ref.  1",  respectively. 

Figure  3 refers  to  the  frontal  point  (-1,  0)  and  shows  results  (a) 
and  (b).  The  real  parts  are  denoted  by  Re,  the  imaginary  parts  by  Im. 
The  phase  angle  is  tt  / 2 - <|>-  While  there  is  sufficient  similarity  between 
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f 


f 


the  two  sets  of  results  (calculated  by  quite  different  methods)  to  create 
confidence  in  the  correctness  of  the  evaluation,  there  are  also  consider- 
able errors  in  (b),in  particular  when  k is  large.  The  same  comment 
applies  to  Fig.  4 where  the  corresponding  results  are  shown  for  the 
collocation  point  on  Eg.  (The  two  curves  for  4>  almost  coincide  in  the 
case  of  Fig.  4 and  are  not  shown). 

(Note  that  the  antisymmetry  of  the  errors  in  Fig.  1 is  not  re- 
peated by  the  unsteady  parts  of  the  errors:  those  in  Fig.  4 are  much 
the  larger.  ) 

Superficially,  one  might  be  tempted  to  conclude  from  Figs.  3 and 
4 that,  by  making  k small  enough  (i.  e.  , by  using  enough  elements  En) 
one  could  make  all  errors  due  to  the  parabolic  approximation  small 
enough  to  make  them  negligible.  However,  this  would  be  a rather  doubt - 

1 

ful  conclusion.  Imagine  a wing  having  a fixed  reference  reduced  frequency 
kc/2‘  "^s  one  decreases  k>  one  decreases  all  the  unsteady  parts  wg^  of 

the  influence  coefficients  themselves;  nevertheless,  the  solution  (the 
unsteady  pressure  distribution)  must  be  almost  independent  of  kand  must 
converge  to  the  correct  solution  as  k -*0.  In  other  words,  as  the  w ^ 

1 

decrease,  their  number  increases.  The  process  k-*-0  is  not  simple  and 
requires  some  detailed  study,  but  it  seems  clear  enough  that  percentage 
errors  should  be  more  directly  instructive  than  the  absolute  errors  that 
are  shown  in  Figs.  3 and  4. 

Percentage  errors  in  the  unsteady  parts  w ^ are  shown  in 
Fig.  5.  For  Re  (x  = +1)  no  percentage  curve  is  drawn  because  here 
there  is  a sign  reversal  in  the  correct  curve.  Fig.  4;  in  this  case,  the 


They  are  usually  entered  as  positive,  regardless  of  their  actual 
sign  (which  can  be  read  from  Figs.  3 and  4).  4 
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error  due  to  the  parabolic  approximation  is  48%  when  k=0.  1,  62%  when 
k=0.2.  The  percentages  that  constitute  the  error  curves  in  Fig.  5 are 
smaller  than  these  values  but  are  still  fairly  large  even  when  k-*-0. 

Also  shown  in  Fig.  5 are  some  "Ref.  1"  results.  The  additional 
error  that  arises  from  the  inaccuracy  in  the  K -values  is  relatively 
small,  particularly  when  k is  large,  but  occasionally  it  can  be  quite 
large,  see  the  point  for  k=0.  1 of  Re*u'(x  = -1).  This  particular  point  is  a 
consequence  of  the  "corner  error"  which  is  discussed  in  Ref.  3.  This  large 
error  is  disturbing  because  k=0.  1 is  well  within  the  range  of  normal  usage. 


As  an  additional  numerical  example  we  compared  for  k=0.  2 
the  unsteady  contributions  to  the  sets  (a)  and  (b)  at  larger  (positive  and 
negative)  values  of  x (up  to  p = 12).  We  found  that  the  errors  (in  (b)  ) 
were  not  large  (they  are  rather  small,  for  k = 0.  2,  even  when  p=l,  see 
Figs.  3 and  4)  but  they  did  not  exhibit  either  the  rapid  convergence  to 
zero  that  is  shown  for  the  steady  case  in  Figs.  1 and  2.  Rather,  in  the 
unsteady  example  the  percentage  errors  tended  to  stay  about  constant. 
Thus,  while  in  the  steady  case  the  parabolic  approximation  was  seen  to 
be  quite  adequate  for  all  collocation  points  sufficiently  removed  from  the 
bound  vortex,  this  desirable  property  is  less  pronounced,  at  best,  in 
the  unsteady  case. 

THE  CASE  y / 0 


If  y/0,  r\  in  the  denominator  of  the  inner  integral  in  Eq.  (9)  has 


to  be  replaced  by  (q-y)  , and 


replaces 


1 


. dq 


0 
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It  is  then  easily  shown  that  (read  y as  |y|  ' 

Ws(x,y;k)  = <y£r;k(y+l))-  -^lw8  <Fl;k(y‘ 1})  ] <18> 

The  case  y/0  is  thus  formally  reduced  to  the  already  resolved  case  y=3. 

The  fact  that  two  coefficients  (on  the  right)  are  required  to  form 
one  (on  the  left)  is  not  a duplication  of  effort,  since  those  on  the  right 
will  each  be  required  with  two  different  values  of  y.  Of  more  concern  is 
the  fact  that  the  effective  reduced  frequencies  on  the  right,  e.g.  k(y+l), 
may  become  larger  (except  if  the  wing  aspect  ratio  is  small)  than  the  limit 
k^  1 for  which  our  Tables  1 and  2 have  been  prepared.  However,  these 
sequences  of  coefficients  exhibit  rapidly  increasing  convergence  as  they 
are  lengthened;  we  found  that  going  to  k = 3 required  only  about  two  more 
terms.  F urthermore,  the  high  accuracy  provided  may  be  more  than  is 
required* 

In  many  cases,  the  present  tables  will  be  sufficient.  Extension  of 
the  tables,  if  required,  is  trivial.  At  worst,  one  can  make  more  extensive 
use  of  direct  numerical  integrations. 

REMARK  ON  EXTRAPOLATION 

The  overall  goal  of  the  present  investigation  (an  earlier  result  of 
which  was  Ref.  3)  is  to  calculate  reliable  reference  solutions  for  certain 
problems  in  harmonic  lifting  surface  theory.  We  aim  to  use  finite 
elements  and  to  perform  extrapolations  k — 0;  this  approach  was  provoked 
by  unexpectedly  favorable  results  in  steady  flow  analyses.  As  an  example, 
consider  a square  wing,  divided  into  N equal  square  elements  E . Use 


15 


two  different  types  of  sending  elements  Eg,  the  Falkner  horseshoe  vortex 
("F")  and  the  Woodward^  element  ("W")  (here  the  lift  is  evenly  distributed 
over  Eg;  however,  we  used  collocation  points  at  85%  chord  while  Woodward 
had  proposed  95%).  With  each  Eg,  solve  the  (steady)  problem  of  the  flat 
square  wing  at  incidence  for  three  values  N and  extrapolate  as  follows: 


CL  04,  N -*■  oo  ~ CLa  ~ CLO,  N 


+ (a/N)+  (b/lO 


In  each  case,  the  three  results  determine  the  three  unknown  constants 
Cj^,  a and  b.  We  found 


CLa  = 1.4597  ("F")  and  1.4602  ("W")  with  N = 4,  6 and  8 
and  found 

CT  „ = 1. 4600  ("F")  and  1. 4601  ("W")  with  N = 6,  8 and  10. 

J_j  Ci 

This  is  a remarkable  success.  In  particular  concerning  the  element 
"FM,  numerous  investigations  have  been  made  relating  to  2-D  flow 
(where  exact  solutions  are  available  for  comparisons).  Here  we  see  a 
result  for  a small  aspect  ratio  wing,  a result  obtained  with  little  numeri- 
cal effort,  which  shows  that  the  two  rather  different  types  of  elements 
Eg  lead  to  the  same  result  (the  present  result  agrees  with,  but  is  more 
accurate  than,  any  of  the  lift  slope  values  for  the  square  wing  found  in 
the  literature).  This  observation  serves  to  establish  the  validity  of  the 
finite  element  approach  in  lifting  surface  analyses. 

The  present  investigation  indicates  that  the  same  extrapolative 
dp^roach  may  involve  complications  in  the  case  of  unsteady  flow.  Since 
one  will  keep  the  reference  frequency  kcy  ^ constant,  the  element  frequency 
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k will  decrease  like  l/  N as  N is  increased,  and  the  unsteady  parts  w gU 

of  the  coefficients  wg  will  decrease  correspondingly  while  the  steady  parts 
( s ) 

w remain  constant.  This  occurs  in  fz~e  of  the  goal  to  determine  the 
s 


specific  effects  of  the  wc 


Furthermore,  the  log  k term  in  Ci(k)  will 


affect  the  limit  behavior  as  N -*■  oo. 

It  thus  appears  that  the  goal  oi  calculating  reliable  solutions  to 
problems  in  unsteady  lifting  surface  theory  requires  further  analytical 
considerations  and  numerical  experimentation. 


CONCLUSION 

To  solve  problems  of  harmonic  lifting  surface  theory  by  means 
of  finite  element  methods  requires  the  determination  of  influence  co- 
efficients and,  in  particular,  requires  evaluation  of  singular  integrals. 
An  accurate  method  of  performing  this  task  has  been  demonstrated. 
Current  methods  (like  the  "parabolic  approximation"  to  the  kernel  K of 
the  integral)  may  lead  to  sizable  errors.  Further  errors  which  arise 
from  inaccurately  evaluating  K are  usually  relatively  small  but  may  be- 
come excessive  when  the  reduced  frequency  k is  small. 
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Table  1.  Coefficients  in  Remainder  Terms,  Eqs.  (14b,  c) 


pTi 1 

P = 2 

mmm 

p > 3 

>1 

o 

0. 042893 

0.  003483 

0. 000732 

l/l6p4  - l/32p6  + 5/256p8 

*1 

0. 001227 

0.  000092 

0. 000019 

1/640  p4 

X2 

0.  000020 

0.  000001 

0 

0 

Bo 

0. 016420 

0.  002428 

0.000747 

l/48p3  - 1/160P5 

B1 

0. 000448 

0.  000062 

0.000019 

1/I920p3 

B, 

0. 000007 

0.  000001 

0 

0 

C, 

Table  2.  Elementary  Integrals  1^  t ggj  (17a) 

I0  = ft  = +1.414214 

^ = |log(l  + V7)  = +0.440687 

= (2Ij-  l/2)/48  = -0.011101 

13  = (61  x - /2)/  5760  = +0.000214 

14  = (30I1-13V'2)/1935360  = -0.000003 
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